Vapor-liquid surface tension of strong short-range Yukawa fluid 
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The thermodynamic properties of strong short-range attractive Yukawa fluids, k = 10, 9, 8 and 7, 
are determined by combining the slab technique with the standard and the replica exchange Monte 
Carlo (REMC) methods. A good agreement was found among the coexistence curves of these 
systems calculated by REMC and those previously reported in the literature. However, REMC 
allows exploring the coexistence at lower temperatures, where dynamics turns glassy. To obtain 
the surface tension we employed, for both methods, a procedure that yields the pressure tensor 
components for discontinuous potentials. The surface tension results obtained by the standard MC 
and REMC techniques are in good agreement. 



I. INTRODUCTION 

The short-range hard-core attractive Yukawa (HCAY) 
fluid has been widely used as a simple model for test- 
ing a variety of theories^— and as a reference system for 
modeling the behavior of real fluids, such as colloidal 
suspensions 7 ^ and protein solutions*^—. Its high popu- 
larity is due to the fact that it captures the most impor- 
tant features of these systems, such as coexistence and in- 
terfacial properties which rule many industrial processes 
and biological phenomena. The HCAY pair potential is 
given by 



u(r) 



, cxp[— k{t — a) /a] 
~ r/a ' 



for r < a, 
for a < r, 



(1) 



where k is the interaction range parameter, e is its depth, 
and a = 1 is the hard core diameter. 

The vapor-liquid phase diagrams of the HCAY flu- 
ids have been accessed by using different simulation 
techniques as well as different theoretical approaches. 
In recent previous works the coexistence properties of 
the HCAY fluid were reported for a wide interaction 
ranged—. For long range interactions, i. e. for small 
k, the reported values obtained by simulation and differ- 
ent theoretical approaches of the coexistence properties 
agree. However, for short-range attractions, there ap- 
pear strong differences between theory and simulations. 
Moreover, not even different simulation approaches lead 
to a good agreement. On the other hand, short-range 
attractive potentials are specially relevant for modeling 
protein solutions, i. e., systems having a key role in bio- 
logical application !] 17 ' 18 . Hence, there is an increasing 
attention on this type of interactions^—. For these 



type of interactions, the vapor-liquid coexistence curve 
turns mctastable since it locates within the vapor-solid 
coexistence curv o 11 ' 25 " — . The simulation of this type of 
systems is rather difficult for standard simulation tech- 
niques. 

Determining the surface tension is by far more com- 
putationally demanding than obtaining the coexistence 
curves. This turns even worse when the potential is 
extremely short-range and discontinuous. This explains 
why the vapor-liquid interfacial properties of discontinu- 
ous potentials have been rarely studied so far. For that 
reason there is growing interest in developing new sim- 
ulation techniques for the efficient computation of the 
surface tension. In this direction, some researchers im- 
plemented clever simulation approaches 2 ^ - — . These ap- 
proaches, however, provide a considerable variation over 
their surface tension results^. Furthermore, there are 
very little data when the potential becomes strong short- 
range. A recent work of SingbJI deals with this issue. He 
reported the interfacial properties of strong short-range 
HCAY fluid with k = 8, 9, 10, using Grand Canonical 
transition- matrix Monte Carlo (GC-TMMC). Anyway, 
the surface tension data are scarce and restricted to the 
temperature region close to the critical point. This lack 
of data at lower temperatures motivated us to implement 
REMC, which should help sampling and thus providing 
new data in this region. 

In view of the previous paragraphs, we understand that 
there is a demand for the improvement of the simula- 
tion approaches, specially for short-range interaction po- 
tentials and at low temperatures. Hence, we combine 
the slab technique with the REMC method to study the 
HCAY fluid for k = 10, 9, 8 and 7, and for a wide temper- 
ature range. The main purposes of this work are three: 
First, to show that the methodology that has been pro- 
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posed to calculate the surface tension for discontinuous 
potential o 29 ' 38 is valid even for strong short-range sys- 
tems. Second, to show that the vapor-liquid coexistence 
and interfacial properties can be obtained at lower tem- 
peratures by using the REMC technique. Third, to re- 
port new thermodynamic properties of the strong short- 
range HCAY fluid for a broader range of thermodynamic 
conditions, and to compare the results with those previ- 
ously reported in the literature when available. 

II. METHODS 

In the following subsection we summarizes the interface 
tension evaluation by the virial route for discontinuous 
potentials. In the next one, some details are given on the 
implementation of the REMC technique. 

A. Interface tension calculation 

The methodology employed to calculate the interface 
tension of discontinuous potentials is given in detail in 
previous paper a 38 ' 39 . Its main idea is to obtain the pres- 
sure tensor through the derivatives of the potential given 
by [U The derivative of the discontinuous contribution in 
Eq. [T] is given by 

*p. = - k T6(r-*l (2) 
dr 

where 5(r — a) is a <5-function. The evaluation of this 
expression can be performed during a MC simulation 
through 

r/ , Q(x)-@(x-Aa) . „ , , 

S(x) = — — -r -, as Act -s- 0, (3) 

Act 

where O(ar) is the unit step function: 0(x) = for x < 
and 0(x) = 1 for x > 0. In this work the parameter A is 
given the following values: A = 0.005, 0.010, 0.015, 0.020 
and 0.025. 

For the virial route, the surface tension is calculated 
by 

7=y{< P ->-^[( P -> + ( P w>]} ! (4) 

where L z is the length of the simulation box in perpen- 
dicular direction to the interfaces. The factor 1/2 is due 
to the existence of two interfaces in the system. 

B. Replica exchange method and simulation details 

As mentioned, we combine the replica exchange Monte 
Carlo method^Sr— with the slab technique^. The gen- 
eral idea is to simulate M replicas of the original system 
of interest, being each replica at a different temperature, 



so that the exchange of microstates among the ensemble 
cells is allowed (swap moves). In this way, replicas at 
high temperatures travel long distances in phase space, 
whereas low temperature systems perform precise sam- 
pling at local regions of phase space. Hence, by introduc- 
ing these swap trials, a particular replica travels through 
many temperatures allowing it to overcome free energy 
barriers. Particular ensembles are not disturbed but en- 
riched by the different contributions of the M replicas. 
This technique allows to easily explore the coexistence 
regions at low temperature, where other methods freeze. 
The formal justification of the swap trials implies the 
definition of an extended ensemble as follows 

M 

Qcxtcndcd = \\ QnVTi , (5) 
i=l 

where Q NVTt is the partition function of the Canonical 
ensemble of the system at temperature TJ, volume V, and 
particle number N . M is the considered number of repli- 
cas of the system, which equals the number of different 
temperatures defining the extended ensemble. Q extended 
may be sampled by M independent standard NVT-MC 
simulations. However, swap trials can be now introduced 
between the replicas which improves the sampling of the 
low temperature ensembles. The acceptation probability 
for this swap trials (performed between adjacent replicas) 
is given by 

P acc =mm(l,exp[(/^ - - Uj)]) (6) 

where Ui — Uj is the potential energy difference between 
replicas i and j, and — f3j is the difference between the 
reciprocal temperatures i and j. Adjacent temperatures 
should be close enough to provide large exchange accep- 
tance rates between neighboring ensembles. In order to 
take good advantage of the method, the ensemble at the 
larger temperature must assure large jumps in configura- 
tion space, so that the small temperature ensembles can 
sample from disjoint configurations. 

We employed parallelepiped boxes with sides L x = 
L y = 8.0cr and L z = 40. Oct for simulating the vapor-liquid 
interface of the HCAY fluid. Each of them is initially set 
with all particles (N = 1000) randomly placed within the 
slab which is initially surrounded by vacuum^. The cen- 
ter of mass is placed at the box center. Periodic bound- 
ary conditions are set in the three directions. Verlet lists 
are implemented to improve performance. Simulations 
were carried out in the vapor-liquid metastable region, 
so the highest temperature is set below the critical tem- 
perature. Other temperatures are fixed by following a 
geometrically decreasing trend. The replicas are equili- 
brated by 10 7 MC steps, while maximum displacements 
are varied to yield acceptance rates close to 30 %. Long 
displacements trials are also considered. These displace- 
ments are important since they allow large jumps in the 
vapor phase with relatively large acceptance rates while 
naturally performing particle transference trials between 
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FIG. 1. Density profiles, p*{z), of HCAY fluids with k = 10 
for different temperatures by means of REMC. Larger tem- 
peratures produce larger vapor densities and lower liquid den- 
sities. 

both phases. The thermodynamic properties are calcu- 
lated by considering additional 4 x 10 7 configurations. 

The new data reported with standard MC simulations 
were made using the same parameters as in previous 
worksi 3 -"— , i. e., they were performed in a parallelepiped 
cell with sides L z = 40a > L x = L y = 10er and 1500 
particles. It is well known that with these parameters 
finite size effects are avoided^—. The new results were 
obtained with 6 x 10 7 MC steps to get more accurate val- 
ues than in previous work 13 . This new results for k = 7 
are slightly different from those previously reported by 
Duda et. al^. Additional standard MC simulations were 
also obtained with N = 1000 and L x = L y = 8<r for sev- 
eral conditions, leading to the same results. This suggests 
that size effects are less important for strong short-range 
potentials. 

The critical density and temperature were calculated 
by using the rectilinear diameters law and the universal 
value of /3 — 0.325 49 ". Our results are given in dimension- 
less units, as follows: r* = r/a for distance, T* = fc^T/e 
for temperature, p* = pa 3 for density, and 7* = 7<r 2 /e 
for surface tension. 



III. RESULTS AND DISCUSSION 

The vapor-liquid interfacial properties of strong short- 
range attractive Yukawa systems with k = 10,9,8, and 
7 were obtained using two simulation techniques, i. e., 
standard MC and REMC. Following we focus on the ob- 
tained results. 

Fig. [1] shows typical interfacial density profiles for 
strong short-range HCAY fluids with k = 10 obtained 
at different temperatures. Well-defined liquid and va- 
por regions are observed which makes us confident that 
the interfaces are stable. Furthermore, the systems show 
bulk vapor and liquid regions of similar volume. This is 
convenient to generate relatively large bulk regions and 
make sure they are fully developed. As it is expected, a 
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FIG. 2. Phase diagrams of HCAY fluids. Open circles cor- 
respond to REMC whereas squares to standard MC. A good 
agreement is found for all cases. 

temperature increase leads to lower and higher densities 
of liquid and vapor phases, respectively. Besides, one can 
clearly see that the interface width increases with tem- 
perature, pointing out a decrease of the interface tension. 
It should be noted that the REMC method yields very 
smooth density profiles (a good sampling of highly un- 
corrected configurations is performed), which allows to 
obtain precise values of the liquid and vapor coexistence 
densities, p* and p*, by taking average in the different 
regions of the profiles. This is done without the need of 
fitting a hyperbolic function, as commonly used. 

Fig. [2] shows the phase diagrams obtained from REMC 
for the studied systems (open circles). This figure is built 
from the values of p* and p* obtained from the profiles 
at different temperatures and for k = 10, 9, 8 and 7. The 
same figure also includes previously reported data for 
k = 10,9 and 7 (open squares)^. The open squares 
corresponding to k = 8 are new data. All the data rep- 
resented with open squares were obtained by using the 
standard MC method 16 . An excellent agreement is found 
between REMC and standard MC (codes are totally inde- 
pendent from one another). The data also agree well with 
the coexistence curves reported by Singh^ 7 - (as shown in 
ref. 16 ). Note that the REMC technique provides data 
at lower temperature, where the standard methodolo- 
gies have sampling issues as the system dynamics turns 
glassy^. At even lower temperatures the vapor-liquid co- 
existence metastability brakes and crystallization takes 
placed (not shown). For obtaining the critical points 
we employed the exponential fitting (with (3 = 0.325) 44 . 
Critical points also show a good match between standard 
MC and REMC simulations. 

As mentioned in section [Hi the accurate evaluation of 
the surface tension implies determining it as a function 
of the parameter A and performing an extrapolation for 
A — > 0. Fig. [3] shows the surface tension as a function 
of A for k = 10 and different temperatures. In this fig- 
ure one can see the linear behavior of the surface tension 
as a function of A for all studied temperatures. This 
allows to obtain an accurate value of this property by 
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FIG. 3. Surface tension against A for k = 10 and different 
temperatures. Upper curves correspond to lower tempera- 
tures. Lines correspond to linear regressions. 
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FIG. 4. Surface tension of strong short-range HCAY fluid 
as a function of temperature and n. Circles are obtained by 
REMC and squares by the standard MC method. Crosses 
corresponds to data reported by Singh^ 7 -. The inset shows 
the REMC reduced surface tension against the reduced tem- 
perature. 



simply performing a linear extrapolation for A — > 0. It 
should be noted that the slopes of the regressions never 
vanish, meaning that a single A value cannot be used for 
an accurate determination of the surface tension. More- 
over, the slopes turn larger for decreasing temperature 
(see Fig. Thus, determining the surface tension for 
the short-range HCAY potential at low temperatures im- 
plies obtaining it as a function of A and then performing 
a careful extrapolation for A — » 0. On the other hand, 
when comparing the behavior of the surface tension as 
a function of A for the HCAY fluid to the square well 
fluid, we notice that the slopes for the HCAY fluid have 
the opposite sign. This is probably the effect of having a 
single discontinuity instead of the two corresponding to 
the square well potential. Anyway, in both cases the lin- 
ear extrapolation is essential for the accurate evaluation 
of the surface tension of discontinuous potentials by the 
virial route. 

The results of the surface tension as a function of 
the temperature for all considered k values are shown 
in Fig. 0] and given in Table Q] The curves behavior is 



TABLE I. Surface tension values for the strong short-range 
Yukawa fluids. The subscripts are the estimated errors. 
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very similar to that found with longer interaction ranges, 
i. e., the surface tension raises as temperature decreases 
and the curves shift to the left with increasing n. As in 
previous plots, circles correspond to REMC and squares 
to the standard MC method. The obtained agreement 
is remarkable. This provides confidence in the proposed 
method for the surface tension evaluatio n 29 i 38 ' 46 of these 
particular extremely short-range and discontinuous po- 
tentials by the virial method. 

We also included in Fig.|4]the data recently reported by 
Singh by means of GC-TMMC method 37 . A comparison 
to these data reveals a good qualitative agreement. That 
is, trends show a good match. Nonetheless, our data 
are systematically above Singh's values, specially close 
to the critical point where differences are above 30%. 
It should be mentioned that differences are consider- 
ably larger than our error bars, which are always smaller 
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than twice the symbol size for both cases (see Table 1). 
These differences in the surface tension values are despite 
the fact that a very good match was found between the 
vapor- liquid phase diagrams obtained by GC-TMMC and 
the standard MC technique, as it was shown in a previ- 
ous worklS. Our belief is that the GC-TMMC technique 
would probably yield larger surface tension data for a 
considerably larger number of MC steps (the author re- 
ported a short running time calculation, 48 — 72hrs). 
Note that we performed 6 x 10 7 and 4 x 10 7 steps with 
the standard MC and the REMC method, respectively, 
to produce the data shown in Fig. H A M = 12 REMC 
run takes approximately a month in a quad core desktop 
machine. Thus, in both large number of steps 

was required. In that sense, the gain of employing REMC 
instead of a standard MC technique for the surface ten- 
sion determination is not very remarkable (we expected 
a much better performance of REMC). A similar finding 
was obtained by de Miguel 3 -^ when comparing the ex- 
panded ensemble approach (a clever method also based 
on the introduction of an extended ensemble and thus 
related to REMC) with the virial route for obtaining the 
surface tension. 

Previous findings for intermediate and long interac- 
tion range of the same pair potential show the surface 
tension and the coexistence curves yield a master curve 
when reduced using their critical parameters^. Further- 
more, a single master curve was also found for the co- 
existence data of the strong short-range HCAY 16 . How- 
ever, the reduced surface tension, 7* = 7 / '(p* 2 ^ 3 T*), ob- 
tained by REMC as a function of the reduced tempera- 



ture, = T* /T*, shown in the inset of Fig. [4] forms an 
imperfect master curve, where slight deviations are ob- 
served. Thus, we cannot assure the obeying of the cor- 
responding states law of the strong short-range HCAY 
based on this property. 

IV. CONCLUSIONS 

The presented work fulfills the three purposes stated 
in the introduction. That is, it shows that the virial 
route calculation of the surface tension for discontinuous 
potential o 29 ' 38 produces reliable data for strong short- 
range systems. For this purpose, the surface tension must 
be obtained as a function of A and then a careful extrapo- 
lation for A — ¥ must be implemented. Then, it reveals 
the utility of the REMC technique to assist the virial 
route for obtaining data at low temperatures. Finally, it 
reports new interface tension values for the strong short- 
range HCAY fluids which contribute to the knowledge 
on the behavior of these systems for a broader range of 
thermodynamic conditions. These new data may help 
testing the results of the existing and forthcoming theo- 
retical approaches. 
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